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We study the equilibrium and non-equilibrium properties of strongly interacting bosons on a lattice in pres- 
ence of a random bounded disorder potential. Using a Gutzwiller projected variational technique, we study the 
equilibrium phase diagram of the disordered Bose Hubbard model and obtain the Mott insulator, Bose glass and 
superfluid phases. We also study the non equilibrium response of the system under a periodic temporal drive 
where, starting from the superfluid phase, the hopping parameter is ramped down linearly in time, and back to 
its initial value. We study the density of excitations created, the change in the superfluid order parameter and 
the energy pumped into the system in this process as a function of the inverse ramp rate r. For the clean case 
the density of excitations goes to a constant, while the order parameter and energy relaxes as 1 /r and 1/r 2 
respectively. With disorder, the excitation density decays exponentially with r, with the decay rate increasing 
with the disorder, to an asymptotic value independent of the disorder. The energy and change in order parameter 
also decrease as r is increased. 

PACS numbers: 



I. INTRODUCTION 

The model of bosons on a lattice interacting repulsively 
through a local interaction in the background of a random 
one-body disorder potential (or the disordered Bose Hubbard 
model [1]) has been used as a paradigm for superfluid in- 
sulator transition in a host of disordered quantum systems. 
These encompass a number of different condensed matter sys- 
tems, from A He on disordered substrates J2) or in porous 
media [3], to dirty superconducting films [4 J and Josephson 
junction arrays [5], to disordered quantum magnets (6J- In 
fact, this model has often been used to describe the relevant 
bosonic degrees of freedom near phase transitions in strongly 
disordered systems. There are three main ingredients in this 
model: a hopping or kinetic energy term for bosons, which 
tend to favour delocalized superfluid phases, an onsite repul- 
sion which tries to localize the bosons to create a Mott insula- 
tor, and an onsite one-body random disorder potential which 
scatters the bosons and lead to loss of coherence of the super- 
fluid. The interplay of these three different terms produces a 
rich phenomenology in these systems, both in its equilibrium 
properties and in terms of non equilibrium dynamics in these 
systems. 

Beyond the traditionally material based phenomena, for 
which it serves as a paradigm, the disordered Bose Hubbard 
model can be realized with ultracold atomic systems EH9), 
which have emerged as the new platform to study the be- 
haviour of model many -body Hamiltonians used in condensed 
matter physics and elsewhere ifTUl . The easy tunability of im- 
plemented Hamiltonian parameters and almost complete iso- 
lation from external environment makes these systems attrac- 
tive candidates to simulate strongly interacting quantum many 
body Hamiltonians, both on the lattice and in the continuum. 
Although cold-atomic systems on optical lattices are gener- 
ally free of disorder (which is inevitably present in solid state 
systems), disorder can be added in a controlled manner either 
by use of speckle potentials (8]|9l or by the use of multiple op- 



tical lattice beams with incommensurate wavelengths ll7l [TTI . 
In either case, the disorder potential (or its distribution in the 
case of speckle potentials) is well characterized and the pa- 
rameters characterizing the disorder potential can be changed 
in a controlled way, in contrast to condensed matter systems, 
where the disorder parameters are unknown a-priori, and are 
mostly determined through a post-hoc process of matching ex- 
perimentally measured quantities (like transport co-efficients) 
to theoretical model calculations. The possibility of controlled 
addition of disorder, thus, makes cold atoms uniquely suited 
to study the effects of disorder on strongly interacting quan- 
tum many-body systems. 

Cold atom systems also provide an added advantage of 
easy access to the internal nonequilibrium dynamics of iso- 
lated interacting systems. The low energy scales (in the ab- 
solute sense), the easy tunability of the Hamiltonian param- 
eters and the almost complete isolation of the system from 
external environment make it very easy to perturb the sys- 
tem from its equilibrium state in a well characterized way and 
then follow the dynamics of the system without the help of 
ultrafast probes. This has opened up the possibility of study- 
ing the quantum dynamics of these systems out of equilib- 
rium 02 ELD • 

Since the early work of Fisher et al 01, me equilibrium 
properties of the disordered Bose Hubbard model has been 
treated with various levels of sophistication from mean field 
theory |[T4l4T6ll to strong coupling expansions ifTTl to Monte 
Carlo techniques [ 18— 20|. In this paper, we provide an al- 
ternative approach to studying the disordered Bose Hubbard 
model based on variational wavefunctions. Our approach is 
applicable in the strongly interacting limit of the model, but 
does not place any constraint on the strength of disorder po- 
tential. The variational approach uses a canonical transforma- 
tion to systematically eliminate processes connecting states 
with large energy difference (~ U, the onsite Hubbard repul- 
sion, or more) and generates an effective low energy Hamil- 
tonian for the system in the strongly interacting limit. This 
effective Hamiltonian is then treated with a Gutzwiller mean 
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field wavefunction. There are several benefits to this approach 
over other standard approaches : i) It captures the strong cor- 
relations generated by the boson repulsion more accurately 
than mean field theory ii) The requirement of disorder aver- 
age makes the problem numerically very resource intensive to 
treat beyond mean field theory. Our semi-analytic approach 
lessens the numerical burden, while keeping essential "beyond 
mean field" correlations, iii) Since this approach generates an 
effective low energy Hamiltonian, it can be easily modified to 
study quantum dynamics in these systems. This is a crucial 
aspect of this approach, which makes it qualitatively different 
from more sophisticated Monte Carlo techniques, specially in 
larger than one dimensions. 



In this paper, we first study the equilibrium phase diagram 
of the 2D disordered Bose Hubbard model on a square lat- 
tice within our approach as a function of Uj J and the chemi- 
cal potential [i for different values of the disorder strength V. 
This yields three phases: (a) an incompressible phase incoher- 
ent Mott insulating phase at large interaction strength, whose 
area decreases with increasing disorder strength (b) a super- 
fluid phase, with coherent condensation of the bosons into a 
single quantum state at small interaction strength, and (c) a 
Bose glass phase in between them, where the system is com- 
pressible, but the phase coherence of the bosonic condensate 
is completely destroyed. We also study the non-equilibrium 
dynamics of the system under the following conditions: the 
system is initialized in its ground state in the superfluid phase. 
The interaction parameter is ramped up linearly in time to a 
very high value and then ramped back linearly to its initial 
value. At the end of this process, we study the density of 
excitations produced in the system, the energy pumped into 
the system and the deviation of the superfluid order parame- 
ter from its initial value, as a function of the rate of the ramp, 
1 jr. In the clean case, the excitation density goes to a con- 
stant, while the order parameter deviation and energy scales 
as 1/t and 1/t 2 in the large r limit. With disorder, the ex- 
citation density shows an exponential decay. The energy and 
order parameter deviation also decreases with increasing r, 
although a scaling form is hard to obtain due to inherent noise 
in the data. 



The paper is organized as follows: In section |TTJ we 
present our variational wavefunction approach and introduce 



the canonical transformation. Section III presents the details 
of obtaining the canonical transformation operator and the ef- 



fective low energy Hamiltonian. In section IV we present the 
equilibrium phase diagram calculated within our approach. In 
section [V]we present the results for the non equilibrium dy- 
namics in the system. Finally, we conclude in section VI with 
a summary of our results and a discussion of limitations of the 
present formalism and ways to improve it. 



II. VARIATIONAL WAVEFUNCTION 

The Hamiltonian of the disordered Bose Hubbard model on 
a square lattice is given by 

H = -Jj2 b l b i + 2 H - X ) + £( w < _ W 

(ij) i i 

where b\ creates a boson on site i and n is the boson num- 
ber operator. Here J is the nearest neighbour hopping energy 
scale, U the on-site Hubbard repulsion, fi the chemical po- 
tential and Vi is the random local potential on site i. Vi is a 
spatially uncorrected random variable drawn from a uniform 
distribution in the range — V/2 < vi < V/2, where V sets the 
energy scale for disorder effects. 

The clean Bose Hubbard model (vi = 0) has a quantum 
phase transition between a strongly interacting incompress- 
ible Mott insulating phase with commensurate integer filling 
at small J/U and a phase coherent superfluid state with large 
number fluctuations at large J/U. The transition is character- 
ized by the vanishing of both the superfluid stiffness and the 
compressibility as one reaches the Mott phase. In the presence 
of disorder, there is an intervening Bose glass phase where the 
superfluid stiffness vanishes, but the compressibility remains 
finite. 

We wish to study the equilibrium phases and dynamics 
in the disordered Bose Hubbard model through a variational 
wavefunction approach. For the equilibrium phase diagram at 
T = 0, we use a variational ground state wavefunction of the 
form 



(2) 



Here \ip ) is a Gutzwiller type local mean field state with vari- 
ational parameters f n i, which satisfies J2 n \ fni\ 2 = 1 to en- 
sure normalization of the state, \n)i is the number state with 
n bosons on site i, and e~ lS is a canonical transformation that 
builds in non-local correlations in the proximity of a Mott in- 
sulator. 

The canonical transformation approach has a long history 
of use in the context of Fermi Hubbard model in the strongly 
interacting limit, where it is used to convert the Hubbard 
model to the so called "t-J" model used in the study of high 
temperature superconductors EH . Recently this approach has 
been adapted successfully to study the equilibrium phases of 
and quantum dynamics in clean Bose Hubbard model Il22ll23l . 
The canonical transform uses the local number states (which 
are eigenstates of the local part of the Hamiltonian) as the 
starting point. Note that in the present formulation of the 
canonical transformation, we do not use a particular local state 
as our starting point (except assuming a local number state), as 
is done in Ref. [22] where the atomic limit Mott phase ground 
state with the same number of particles on each site is used 
as the starting point. The hopping terms then start to build in 
correlations between different number states on neighbouring 
sites. 

The hopping term can connect local number states which 
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differ in energy by ~ U or higher. To see this consider a state 
with m particles on site i and w,2 particles on site j and a 
hopping process where a particle hops from j to i. The energy 
difference (coming from the local part of the Hamiltonian) be- 
tween the initial and final state is Se = U(ni—n,2+l)+Vi—Vj 
and \Se\ can be ~ U or more depending on n\ and rz.2. We 
would like to note that, since we are interested in energy 
difference of states connected by hopping (which does not 
change the total number of particles in the system), the chem- 
ical potential drops out of the expression for the canonical 
transform. Hence our formalism is applicable for any [i, even 
to the parameter regime where /i ~ nU, n being an integer. 

The basic idea of the canonical transform is to eliminate 
terms in the Hamiltonian which connects local number states 
differing by a large energy (~ U or more) order by order in 
J/U through the canonical transformation. The easiest way 
to see this is to note that for any operator A, 



(tl>\A\fl>)=(i/) \A*\fo), where A* 



Js Ae -iS 



(3) 



is the canonically transformed operator. For the Hamiltonian, 
H, the requirement that H* does not have any terms con- 
necting states which differ in energy by ~ U fixes the form 
for S. The low energy effective Hamiltonian, H* , obtained 
by the canonical transform, not only allows the low energy 
hopping processes, but also builds in correlations from vir- 
tual transitions to high energy states. In the next section, we 
provide the details of the canonical transformation and the ef- 
fective Hamiltonian for the disordered Bose Hubbard model. 
We would like to note here that although we will focus here on 
a random disorder potential, our formalism is capable of han- 
dling any one-body potential, e.g.it can be used to treat effects 
of harmonic traps in ultracold atomic gases in optical lattice. 



III. THE CANONICAL TRANSFORMATION 

The disordered Hubbard model can be separated into a local 
part containing the interaction and the one body potential and 
a kinetic energy part. 



u 



- y ^2,n i (n i - V) - fiihi (4) 



where jjLi = fj, — /x being the chemical potential and i>i 
the random disorder potential. It is also easy to see that the 
hopping term = — Jo[bj connects local states differing in 
energy by efj — all + Vi — Vj, where a — 0, ±1, ±2... This 
suggests breaking up the hopping term, = ^ Q Tfj, where 

T?- = -Jb\bjS{ni-nj-a + i) (5) 
= -J ^9a\ n + l )i\ n - a )j j( n -a + l\ 



where g n a = y/(n+l)(n -a + l). Here T£ connects states 
with energy difference . A mathematical way of represent- 



ing this information is the identity 

l U 0, J-iA — Eijl-ij, 



(6) 



which will be useful later in deriving the canonical transfor- 
mation operator. 

For weak disorder (V -C U), it is evident that repre- 
sents a low energy hopping process, while for a ^ 
changes energy of the state by an amount ~ all and has to 
be eliminated by the canonical transform. This breakup of the 
kinetic energy term follows the method of Girvin et al ||2~T1 for 
fermionic Hubbard model with one crucial difference: in the 
Fermi Hubbard model, the local Hilbert space is constrained 
by Pauli exclusion and hence a = 0, ±1, whereas in the 
Bosonic model, the infinite Hilbert space leads to a taking 
all possible integer values. In practice, the local Hilbert space 
is cutoff at some high value of occupancy number, and a will 
be restricted accordingly. This formalism can be generalized 
to strong disorder potentials with some more complications, 
which will be discussed in a future work. 

The canonical transformation operator iS has an expansion 
in J/U, i.e. iS = iS 1 + iS 2 + where iS m - (J/U) m 
and terms upto iS m completely removes high energy terms 
upto order J( J/U)" 1 ^ 1 . Using the identity, eqn. (|6h, it can be 
shown that 



^ - EE 



(7) 



(ij) y 

removes all high energy terms O(J), while 



(8) 



i E E 



( F -P , F a \ 
(ij)(kl) a^P^O \ b kl ' ^u) 



ij) l £ ij £ kl 



removes high energy terms upto 0(J 2 /U). 
The effective Hamiltonian H* is then given by 



(9) 



{ij) 



(ij)(kl) a^O 



which, in the clean case of a Bose Hubbard model without 
disorder, reduces to 



H*(V = 0)=H 



E E 

(ij) (kl) Q#0 



rp-l 

ij' 1 hi 



aU 



(10) 



The effective low energy Hamiltonian thus consists of three 
terms : (a) Hq which gives the local interaction and disorder 
potential, (b) T , which represents the low energy hopping, 
and (c) the last commutator, which can be easily interpreted 
as a second order perturbation, and takes care of virtual tran- 
sitions to high energy states. 

We would like to note that our canonical transformation im- 
proves upon previous formulation by Trefzger et. al EH l23l 
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FIG. 1: The zero temperature equilibrium phase diagram of the 
disordered Bose Hubbard model in the n/U-J/U plane for (a) 
V/U — 0.3 and (b) V/U = 0.6 respectively. The phase to the right 
of the thick red line is the superfluid phase, while the Mott phase is 
enclosed by the dotted blue line. The phase in between is the Bose 
glass phase. 



in the following ways: i) It can handle non-uniform states with 
arbitrary one-body potentials in the local part of the Hamilto- 
nian, which is crucial in treating disordered bosons and ii) It 
takes into account the full Hilbert space for bosons and is not 
an expansion around a state with a fixed number of particles 
on each site. This is crucial to look at the Bose glass phase 
(and to study properties of bosons in a trap), where the den- 
sity varies from site to site. This formulation can also handle 



more accurately the superfluid phase near the Mott lobes, as 
it treats all the states in the local Hilbert space on equal foot- 
ing. In fact, in the clean case, if one keeps only three states in 
the local Hilbert space (the commensurate density in the Mott 
state, no an d n o ± 1), then a is restricted to 0, ±1, and our 
formulation reduces to that in Ref . [22] 



IV. EQUILIBRIUM PHASE DIAGRAM 

Starting from the early prediction of Fisher et. al |[T), 
the equilibrium phase diagram of disordered Bose Hubbard 
model has been worked out by several previous authors us- 
ing various techniques ranging from mean field theory [14] to 
quantum Monte Carlo techniques [18 19] . Although our main 
motivation is to study dynamics of the system when interac- 
tion parameters are tuned, we present the equilibrium phase 
diagram obtained by our method for the sake of completeness. 
This will also set the stage for our study of dynamics in two 
ways: (i) the equilibrium ground state forms the initial condi- 
tion for the dynamics of the system and (ii) we would like to 
know the trajectory of the system, i.e. whether it goes into the 
Mott or the Bose glass phase as we ramp up the interaction 
parameter starting from the superfluid phase. 

The ground state is obtained by minimizing the energy in 
the variational state, which is equivalent to minimizing the 
expectation of H* in the the mean field state |^o)- A straight- 
forward algebra shows that the ground state energy is a sum 



of six different contributions, £ = J2 r =o ^"»"> wnere 



— n(n 
2 v 



(ii) 



is the local energy corresponding to the interaction and disor- 
der potential, 



£l = - J XI ( n + Vfn+lifnifnjfn+lj 
n{ij) 

is the low energy nearest neighbor hopping, 



(12) 



£ 2 _ — ^ J2(n+l)\f |2 ( n + a + l)l/n+q+ij| 2 Z (»-« + l)l/»-q+ij| 2 



(13) 



(ij) na^O 

is the second order density-density interaction energy, 

J 



n+2j 



(ij) 



1 1 



(14) 



is a second order pair hopping term where two bosons hop to the nearest neighbor. 

j2 f* f n . 

£4 = y X £ «+2» m [ g n a g2+if* +ak f n+a+lk f*_ a+lj f n _ a+2j - [a _> -a)] + h.C. 



(15) 
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which represents a second order process where two bosons from two different neighboring sites hop onto a site and its reverse 
process, and finally 

(ij)(jk) na^Q a 

which represents a second order next nearest neighbor hopping process. 

I 



The energy is then minimized with respect to the varia- 
tional parameters f ni to obtain the ground state wavefunction. 
The three different phases are then identified according to the 
following criterion: The superfiuid phase is characterized by 
a non vanishing superfiuid stiffness, which controls the en- 
ergy of the system for long wavelength distortion of the phase 
of the Bose Einstein condensate. This can alternatively be 
thought as the diamagnetic response of the system to a vector 
potential. In presence of a static vector potential A along the x 
direction, the hopping parameters acquire an Ahronov-Bohm 
phase, Jij — > Jije A ( Xi ~ Xi \ and correspondingly H — > Ha 
and S — > Sa- The superfiuid stiffness can then be calculated 
as 

_ 1 v d 2 TOo 

P'-X^~d^~ lA -° (17) 

c 

where C denotes disorder configurations and N c is the number 
of configurations kept in the disorder average (typically ~ 100 
in our calculations). We note that in case of finite disorder 
potential we will always work with disorder averaged quan- 
tities in this paper. Any state with p s ^ will be identified 
as a superfiuid phase. In the non-superfluid phase, we distin- 
guish between the Bose glass phase and the Mott insulating 
phase by the fact that the Bose glass phase has a finite com- 
pressibility, while the Mott insulating phase is incompressible. 
Within our formalism, this implies that /„ j = 1 for all the 
lattice sites in all the disorder configurations in commensurate 
Mott insulator of filling n , while in the Bose glass phase, 
max(f no i) < 1. We note here that although f no i = 1 for all 
lattice sites in a Mott phase, the canonical transform mixes in 
virtual number fluctuations in the ground state wavefunction. 

In Fig. [T] (a) and (b), we study the phase diagram of the 
system in the J/U-fi/U plane, focusing in and around the 
no = 1 Mott plateau, for different disorder strengths V/U — 
0.3 and V/U — 0.6 respectively. The parameter regime to the 
right of the thick red line represents the superfiuid phase with 
a non-zero superfiuid stiffness (p s =^ 0). The region enclosed 
to the left of the blue dotted line is the incompressible Mott 
phase, while the region in between these two lines represents 
the Bose glass phase with non zero compressibility but zero 
superfiuid stiffness. 

The phase diagram qualitatively captures the basic physics 
of the disordered Hubbard model. In the atomic limit, (J = 
0), the system remains in the Mott phase as long as V/2 < 
fi < U — V/2. The local Hamiltonian Ho is then optimized 
by the configuration of one particle on each site. On the other 
hand, for /i < V/2, there are sites where the local Hamilto- 
nian is optimized by a hole, while for /j, > U — V/2, there are 



sites where the local Hamiltonian is optimized by double oc- 
cupancy. Thus the state in this limit has number fluctuations 
(and hence is compressible) while the local nature of the fluc- 
tuations imply that superfiuid stiffness is 0. This state is thus 
in the Bose glass phase. As J/U is increased, the Mott phase 
first gives rise to a narrow region of Bose glass phase, which 
then gives way to the superfiuid phase. In the region, where 
he atomic limit ground state is a Bose glass, we see a direct 
transition between a Bose glass and a superfiuid phase. 

With increasing disorder strength, we find two distinct fea- 
ture of the phase diagram: (a) The Mott region shrinks with in- 
creasing V/U and (b) The direct Bose glass to superfiuid tran- 
sition takes place at larger values of J/U. The first one can 
be easily explained by noting that with increasing V/U, the 
width of the Mott phase in the atomic limit (U — V/2 > /j, > 
V/2) decreases. The second feature is explained by the fact 
that stronger disorder leads to stronger scattering and hence 
larger values of J/U is required to restore phase coherence 
and hence superfluidity in the system. 

Before concluding this section, we note that analogous 
phase diagrams, which are in qualitative agreement with ours, 
have been derived using single-site and multi-site mean-field 
mean-field theories |[T4l4T6ll . strong-coupling expansions ifTTl 
and quantum Monte Carlo lfT8H20l . All of these methods con- 
cur with ours regarding the qualitative features of the phase 
diagram such as presence of a glassy region between the su- 
perfiuid and Mott phase in the presence of disorder and the 
increase in the extent of this glassy region towards the edge of 
the Mott lobes. 



V. DYNAMICS IN THE DISORDERED BOSE HUBBARD 
MODEL 

Our main goal in this paper is to study the dynamics of the 
disordered Bose Hubbard model when a Hamiltonian parame- 
ter (in our case the hopping J) is changed in time. Although a 
lot of work has been done on the equilibrium phase diagram of 
the disordered Bose Hubbard model, very little is known about 
the dynamics of this system. In this context it is worth noting 
that the variational wavefunction and the canonical transfor- 
mation is especially well suited to treat the dynamics in this 
system. For a dynamically changing system one can write 
down a variational wavefunction of the form 

|^(t)> = e-^ J ^\Mt)) \Mt)) =IIE/™(*)l n >* 

i n 

(18) 
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FIG. 2: Disordered averaged defect density as a function of rate of 
change of hopping for: (a) top panel: clean case with fi/U = 0.5 
(left) and fi/U = 0.9 (right) (b) middle panel: disorder potential 
V/U = 0.4 with fi/U = 0.5 (left) and fi/U = 0.9 (right) (c) lower 
panel: disorder potential V/U = 0.6 with fi/U = 0.5 (left) and 
H/U = 0.9 (right). 



FIG. 3: Relaxation of superfluid order parameter as a function of 
rate of change of hopping. The deviation of |$(2r)/$(0)| from 1 is 
plotted for: (a) top panel: clean case with \i/U = 0.5 (left) and 
fi/U — 0.9 (right) (b) middle panel: disorder potential V/U = 
0.4 with fx/U = 0.5 (left) and fi/U = 0.9 (right) (c) lower panel: 
disorder potential V/U = 0.6 with fi/U = 0.5 (left) and fi/U = 0.9 
(right). 



where the canonical transformation is evaluated with the in- 
stantaneous value of the parameter J(t). Note that since 
our canonical transform was based on the idea of eliminat- 
ing terms in the Hamiltonian, which connects states differing 
by a large energy ~ U, this would lead to a coarse grained 
dynamics valid for timescales much larger than U . Further, 
since the canonical transform was not an expansion around a 
particular state (like the Mott state), this can faithfully cap- 
ture the evolution of the excitations that are inevitably created 
during time evolution. 

The Schrodinger equation can then be written as 

i\i> ) = {H* - S*)\^ ) (19) 

where S* — e tS Se~ zS and the initial ground state is evolved 
according to this equation. 

At this point it is useful to look at the particular form of 
dynamics we are interested in. We start our system in the 
ground state with an initial value of Ji, which puts it in the 
superfluid phase. We then decrease J linearly to a very small 
value Jf close to the atomic limit with a rate r _1 . We then 
ramp back to our initial value Ji with the same ramp rate. The 



explicit time dependence of the hopping parameter is given by 

J(t) = J l + {Jf-J i ) t - t<T 

= J f + (Ji - Jf) 1 - t>r (20) 

The effective Hamiltonian H* gives rise to energy scales of J, 
U and J 2 /U, while the S term generates scales of A J/(Ut), 
AJJ(t)/U 2 T 2 etc. We assume Ut > 1 (later we will mostly 
be interested in the regime Ut ^ 1), and we will only keep 
the first order term S 1 in the dynamical equations. This leads 
to a notable simplification; since iS 1 oc iS , iS* = iS, i.e. 
there is no Berry phase contribution from rotating the iS term. 
We note that this simplification goes away if we include higher 
order terms in iS. 

We are interested in the excitations created as we ramp 
down to the atomic limit and ramp back up to the initial value 
of J/U. To study this we look at the defect density which, for 
a given disorder configuration is given by 

^( r ) = ^E 1 -!^o(2r)|^(0))| 2 (21) 

i 

where l^oW) = fni(t)\n)i is the local Gutzwiller wave- 



7 



function at time t. Note that since the final and initial values 
of J/U are same, the canonical transformation operator does 
not affect this definition of defect density. For the disordered 
Bose Hubbard model we study the defect density averaged 
over many disorder realizations. The defect density is plotted 
as a function of the time constant r for various values of V/U 
and fx/U in Fig[5] The top panel shows the clean case (V — 0) 
results for (left): fi/U = 0.5, where one passes close to the 
Mott lobe tip as one decreases the hopping J, and (right): 
fi/U = 0.9, which is far away from the Mott tip. The defect 
density shows oscillatory behaviour with t with the large r 
(Ut 1) limit exhibiting an envelope which is constant with 
t. We note that we have taken care to fix the gauge during 
the time evolution and hence the oscillations are not a result 
of the system sampling different gauge configurations in time. 
Rather, the return of the system to its initial state (pd = 0) 
for certain rates of change of hopping has similar origins as 
found in Ref. O where the system was found to return to 
its initial state under the influence of a periodic drive for cer- 
tain drive frequencies. The constant envelope characterizes 
the fact that within the canonical transformation there is a low 
energy state orthogonal to the initial ground state (with the de- 
generacy broken on a scale of ~ J 2 /IT), The limiting constant 
value is ~ 1 away from the Mott tip, where it is easy to create 
excitations and goes to ~ 0.6 near the Mott tip. The middle 
panel shows the defect density as a function of the ramp rate 
for a disordered system characterized by V/U = 0.4, while 
the lower panel shows a system with V/U = 0.6. The oscil- 
latory behaviour persists, but presence of disorder damps the 
oscillations, with the defect density showing an exponential 
decay with the inverse ramp rate. Disorder leads to scattering 
and lifts the degeneracy of the low lying state, thus leading 
to an exponential decay of the defect density. It is also clear 
by comparing the middle and the lower panel figures that as 
V/U increases, the damping timescale becomes smaller. In 
the large r limit, the defect density goes to a constant value, 
which is almost independent of V/U and depends crucially 
on fi/U. For fi/U = 0.5, this value is ~ 0.3, while for 
fi/U = 0.9, this value increases to ~ 0.6. The finite value 
of the defect density in the large r limit is expected, as the 
system starts from superfluid phase with associated gapless 
modes and hence there is no excitation gap to protect defect 
creation in the slow ramping limit. 

We have also studied the evolution of the superfluid order 
parameter 



N. 



(22) 



as the hopping is ramped down and up. To ensure normaliza- 
tion it is easiest to look at the ratio r = $(2r) /$(0) and then 
construct the disorder average of this quantity. The disorder 
averaged r goes to 1 in the large r limit and hence we look 
at \r\ — 1 to determine how the order parameter relaxes to its 
initial value as a function of the ramp rate. This quantity is 
plotted in Fig. [3] for the clean case (top panel) and the disor- 
dered case for V/U — 0.4 (middle panel) and V/U = 0.6 
(lower panel). In the clean case, the quantity |r| — 1 clearly 
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FIG. 4: Residual energy in the system as a function of rate of change 
of hopping for: (a) top panel: clean case with fi/U = 0.5 (left) and 
fi/U — 0.9 (right) (b) bottom panel: disorder potential V/U — 0.4 
with n/U = 0.5 (left) and fi/U = 0.9 (right). 



shows a power-law scaling with an asymptotic 1/r envelope 
on top of oscillations in the large t limit. The relaxation in the 
disordered case is not so simple. Although \r\ — 1 goes down 
with r, we have not been able to clearly extract either a power 
law or an exponential scaling from the large r limit. There is 
a substantial window of r values, where a 1/r power law can 
be defined, but the data seems to deviate from this scaling for 
larger values of r. 

We also study the residual energy pumped into the system 
in the process of ramping down the hopping and returning 
back to the original value. This is important as there is a lack 
of in-situ measurements of temperature in optical lattices and 
the energy pumped into the system is often taken as a bound 
on the amount of heating in the system. The excess energy of 
the system in the final state is given by 

Q = (V(2t)|#|V(2t)) - M0)|ff|V(0)) (23) 

The excess energy of the system as a function of ramp rate 
is shown in Fig. |4] In the clean case, the energy decays with 
r, with a 1 /r 2 envelope in the large r limit. One way to un- 
derstand this is that within a Gross-Pitaevsky description, the 
lowest order dependence of the energy on the order parameter 
is E ~ |$| 2 , and so, a 1/t relaxation of the order param- 
eter leads to a 1/r 2 energy relaxation in the system. In the 
disordered case, although the excess energy decreases with r, 
numerical accuracy of the data forbids a clear extraction of an 
asymptotic limit. 



VI. CONCLUSIONS 

In this paper, we have studied the equilibrium and non- 
equilibrium properties of the disordered Bose Hubbard model 
using a new variational wavefunction approach. Our vari- 
ational wavefunction implements the canonical Schrieffer- 
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Wolf transformation, which has been extensively used for 
strongly interacting Fermions to the case of strongly repulsive 
Bosons in a non-uniform potential background. We have de- 
termined the equilibrium phase diagram of the system, which 
shows the expected Mott insulator, Bose glass and the super- 
fluid phases. Our phase diagram is qualitatively similar to the 
phase diagram obtained by more sophisticated techniques. 

We have also studied the non-equilibrium properties of the 
disordered bosons within this variational approach. We have 
focused on the specific dynamic process, where, starting from 
the system in its ground state in the superfiuid phase, the hop- 
ping parameter J is ramped down linearly with a rate r _1 to 
a value very close to 0. The hopping is then ramped back with 
the same rate to its initial value. We look at the response of the 
system to this non-equilibrium cyclic process by studying the 
density of excitations, the superfiuid order parameter and the 
energy of the excitations in the final state obtained from the 
time evolution of the initial ground state. In the clean system 
(V = 0), we find that all the three quantities show oscilla- 
tory behaviour as a function of the inverse ramp rate r. The 
asymptotic envelope of the defect density is a constant, while 
the order parameter and the energy decays as 1/r and 1/r 2 
respectively in the large r (Ut 3> 1) limit. In the disordered 
system, the oscillations persist, although they are damped by 
disorder scattering. The defect density, as a function of the 
inverse ramp rate, oscillates with an exponentially decaying 
envelope. The decay of the defect density as a function of 
r increases with increasing V/U, while the value of the ex- 
citation density in the large r — > oo limit is independent of 



V/U, but depends on the value of fi/U, i.e. it increases as 
one moves away from the tip of the Mott lobe. The deviation 
of the superfiuid order parameter from its initial value as well 
as the energy pumped into the system decays with decreasing 
ramp rate, but numerical noise prohibits a clear extraction of 
a large r asymptotic scaling. 

The variational wavefunction and the associated canonical 
transform used in this paper provides a new analytic way of 
treating the problem of disordered strongly interacting bosons 
on a lattice. In fact, this technique can be used to study 
any one body potential (including harmonic trap potentials 
relevant to the cold atom experiments). In its current form, 
the variational wavefunctions capture the essential physics of 
the superfiuid, Bose glass and Mott insulator phases in the 
low disorder limit V/U <C 1. Although we have stretched 
the technique to V/U ~ 0.6, the numerical accuracy of the 
method decreases and quantitative match of the phase diagram 
with the Monte Carlo results deteriorates. A more complete 
formulation, which is beyond the scope of this paper, would 
incorporate the fact that for U > \ Vi — Vj | J, the hopping 
on the bond between i and j is completely frozen, while for 
Vi — Vj ~ nU, there is a low energy hopping process which 
changes the number of multiple occupancies in the system. 
Further development along these lines would lead to higher 
numerical accuracy and wider applicability of this new tech- 
nique. 

This work is supported by AFOSR JQI-MURI, ARO- 
DARPA-OLE, and NSF-JQI-PFC. 
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